#pragma once

void init_4F_projector();
void Projection_TL_OpPC(O5_ccccssss OpPC,colorcolorspinspin *Sinv_x0, colorcolorspinspin *Sinv_0x,int r1,int r2,int r3,int r4);
void Projection_TL_OmPC(O5_ccccssss OmPC,colorcolorspinspin *Sinv_x0, colorcolorspinspin *Sinv_0x,int r1,int r2,int r3,int r4);
void Projection_TL_OpPV(O5_ccccssss OpPV,colorcolorspinspin *Sinv_x0, colorcolorspinspin *Sinv_0x,int r1,int r2,int r3,int r4);
void Projection_TL_OmPV(O5_ccccssss OmPV,colorcolorspinspin *Sinv_x0, colorcolorspinspin *Sinv_0x,int r1,int r2,int r3,int r4);
void  Projection_TL_OPC( O5_ccccssss OpmPC, colorcolorspinspin Sinv1,colorcolorspinspin Sinv2,colorcolorspinspin Sinv3,colorcolorspinspin Sinv4, int ipm);
void  Projection_TL_OPV( O5_ccccssss OpmPV,colorcolorspinspin Sinv1,colorcolorspinspin Sinv2,colorcolorspinspin Sinv3,colorcolorspinspin Sinv4,int ipm);
void Projection_TL_2r_OPC( r_O5_ccccssss OpmPC, colorcolorspinspin *Sinv_x0,colorcolorspinspin *Sinv_0x,int ipm);
void Projection_TL_2r_OPV(r_O5_ccccssss OpmPV, colorcolorspinspin *Sinv_x0,colorcolorspinspin *Sinv_0x,int ipm);
void Projection_TL_2r_OpPC(r_O5_ccccssss OpPC,colorcolorspinspin *Sinv_xy, colorcolorspinspin *Sinv_yx);


//////////////////////////////////////////
// Notation from 
// arxiv:help-lat/9902030v1
// (bibbia)
///////////////////////////////////////////

// Green function Gpm: (eq 25 bibbia)
// Gpm_g1g2= ps1 bar(psi2) Opm_g1g2 psi3 bar(psi4)
// with (eq 26)
// Opm(0)_g1g2=1/2 [ bar(psi1)(0) g1 psi2(0) bar(psi3)(0) g2 psi4(0) +- bar(psi1)(0) g1 psi4(0) bar(psi3) (0)g2 psi2(0) ]
// Thus
// Gpm_g1g2(0) = (psi_1 (x)bar{psi}_1(0) ) (psi_3 (x)bar{psi}_3(0) ) ( g1 x g2 ) (psi_2 (0)bar{psi}_2(x) ) (psi_4 (0)bar{psi}_4(x) )
// -> Gpm_g1g2 (0) = S1(x,0) S3(x,0) (g1xg2) S2(0,x) S4(0,x)

// In terms of 2F Non amputed GF (eqs 27-28)
// Gpm_g1g2 (0) = G_r1r3*G_r2r4 -+ G_r1r4*G_r2r3 (change of sign due to the grassmann variables
// If r=0001 or r=1110 -> Gpm_g1g2=GOS*GTW -+ GTW*GOS

// So the propagators to ampute are (eq 31)
// Lambda = S1(x,0) ^-1 S3(x,0)^-1 Gpm_g1g2(0) S2(0,x)^-1 S4(0,x)^-1 (note the order r1r3r2r4)

O55 PpPC, PmPC;
O55 PpPV, PmPV;

O5 SS2r, PP2r, VV2r, AA2r, TT2r;
O5 SP2r, PS2r, VA2r, AV2r, TT52r;

O5 SS, PP, VV, AA, TT;
O5 SP, PS, VA, AV, TT5;

int pm(int ipm)
{
	if (ipm == 0 ) return 1;
	else if (ipm == 1 ) return -1;
	else return 0;
}


void init_4F_projector()
{

      O55 *PpmPC=malloc(sizeof(O55)*(2));
      O55 *PpmPV=malloc(sizeof(O55)*(2));

	// eq 37 & 38 from bibbia
      for (int ipm=0; ipm<2; ipm++)
      {
      PpmPC[ipm][0][0]=1./(64.*Nc*(Nc+pm(ipm)));
      PpmPC[ipm][1][1]=1./(64.*(pow(Nc,2)-1.));
      PpmPC[ipm][1][2]=pm(ipm)/(32.*Nc*(pow(Nc,2)-1.));
      PpmPC[ipm][2][1]=pm(ipm)/(32.*Nc*(pow(Nc,2)-1.));
      PpmPC[ipm][2][2]=1./(16.*(pow(Nc,2)-1.));
      PpmPC[ipm][3][3]=(2.*Nc+pm(ipm))/(32.*Nc*(pow(Nc,2)-1.));
      PpmPC[ipm][3][4]=-pm(ipm)/(32.*Nc*(pow(Nc,2)-1.));
      PpmPC[ipm][4][3]=-pm(ipm)/(32.*Nc*(pow(Nc,2)-1.));
      PpmPC[ipm][4][4]=(2.*Nc-pm(ipm))/(96.*Nc*(pow(Nc,2)-1.));

      PpmPV[ipm][0][0]=-1./(64.*Nc*(Nc+pm(ipm)));
      PpmPV[ipm][1][1]=-1./(64.*(pow(Nc,2)-1.));
      PpmPV[ipm][1][2]=-pm(ipm)/(32.*Nc*(pow(Nc,2)-1.));
      PpmPV[ipm][2][1]=-pm(ipm)/(32.*Nc*(pow(Nc,2)-1.)); 
      PpmPV[ipm][2][2]=-1./(16.*(pow(Nc,2)-1.)); 
      PpmPV[ipm][3][3]=(2.*Nc+pm(ipm))/(32.*Nc*(pow(Nc,2)-1.));
      PpmPV[ipm][3][4]=-pm(ipm)/(32.*Nc*(pow(Nc,2)-1.));
      PpmPV[ipm][4][3]=-pm(ipm)/(32.*Nc*(pow(Nc,2)-1.)); 
      PpmPV[ipm][4][4]=(2.*Nc-pm(ipm))/(96.*Nc*(pow(Nc,2)-1.)); 
      }

      for (int i=0; i<5; i++) for (int j=0; j<5; j++){
		PpPC[i][j]=PpmPC[0][i][j];
		PmPC[i][j]=PpmPC[1][i][j];		
                PpPV[i][j]=PpmPV[0][i][j];
                PmPV[i][j]=PpmPV[1][i][j];
	}

}

void Projection_TL_OpPC(O5_ccccssss OpPC,colorcolorspinspin *Sinv_x0, colorcolorspinspin *Sinv_0x,int r1,int r2,int r3,int r4)
{
       // Lambda = S1(x,0) ^-1 S3(x,0)^-1 Gpm S2(0,x)^-1 S4(0,x)^-1
       //  Returns SS, PP, VV, AA, TT which are 
       // the Lambda functions projected onto each GammaXGamma but without the factor PpmPX
       Projection_TL_OPC(OpPC,Sinv_x0[r1],Sinv_x0[r3],Sinv_0x[r2],Sinv_0x[r4], 0);
}

void Projection_TL_OmPC(O5_ccccssss OmPC,colorcolorspinspin *Sinv_x0, colorcolorspinspin *Sinv_0x,int r1,int r2,int r3,int r4)
{
        Projection_TL_OPC(OmPC,Sinv_x0[r1],Sinv_x0[r3],Sinv_0x[r2],Sinv_0x[r4], 1);
}

void Projection_TL_OpPV(O5_ccccssss OpPV,colorcolorspinspin *Sinv_x0, colorcolorspinspin *Sinv_0x,int r1,int r2,int r3,int r4)
{
        Projection_TL_OPV(OpPV,Sinv_x0[r1],Sinv_x0[r3],Sinv_0x[r2],Sinv_0x[r4], 0);
}

void Projection_TL_OmPV(O5_ccccssss OmPV,colorcolorspinspin *Sinv_x0, colorcolorspinspin *Sinv_0x,int r1,int r2,int r3,int r4)
{
        Projection_TL_OPV(OmPV,Sinv_x0[r1],Sinv_x0[r3],Sinv_0x[r2],Sinv_0x[r4], 1);
}



///////////////////////////////////////////////////
//////////////////////////////////////////////////
//////////////////////////////////////////////////
void  Projection_TL_OPC( O5_ccccssss OpmPC, colorcolorspinspin Sinv1,colorcolorspinspin Sinv2,colorcolorspinspin Sinv3,colorcolorspinspin Sinv4, int ipm)
{

        ccccssss SS_temp; ccccssss PP_temp; ccccssss VV_temp; ccccssss AA_temp; ccccssss TT_temp;
        ccccssss SS_temp2; ccccssss PP_temp2; ccccssss VV_temp2; ccccssss AA_temp2; ccccssss TT_temp2;


        int irho, irhop;
        complex Proj_factor;

            for (int iop=0; iop<5; iop++)
            {
                // Initialize
                SS[iop]=0;
                PP[iop]=0;
                VV[iop]=0;
                AA[iop]=0;
                TT[iop]=0;


                //External indices to be contracted with the projector
                for (int rho=0; rho<4; rho++) for (int rhop=0; rhop<4; rhop++)
                for (int R=0; R<3; R++) for (int Rp=0; Rp<3; Rp++)
                {
                        SS_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;
                        PP_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;
                        VV_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;
                        AA_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;
                        TT_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;


                      // External indices to be contracted with the right propagators   
                      for (int beta=0; beta<4; beta++) for (int delta=0; delta<4; delta++)
                      for (int B=0; B<3; B++) for (int D=0; D<3; D++)
                      {
                        SS_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;
                        PP_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;
                        VV_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;
                        AA_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;
                        TT_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;

                        // Internal indices to contract with left propagators
                        for (int alpha=0; alpha<4; alpha++) for (int gam=0; gam<4; gam++)
                        for (int A=0; A<3; A++) for (int C=0; C<3; C++)
                        {
                                // SS
                                Proj_factor=(base_gamma[0]).entr[rho]*(base_gamma[0]).entr[rhop];
                                irho=(base_gamma[0]).pos[rho]; irhop=(base_gamma[0]).pos[rhop];
                                SS_temp[R][B][Rp][D][rho][beta][rhop][delta]=SS_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPC[iop][A][B][C][D][alpha][beta][gam][delta];

                                // PP
                                Proj_factor=(base_gamma[5]).entr[rho]*(base_gamma[5]).entr[rhop];
                                irho=(base_gamma[5]).pos[rho]; irhop=(base_gamma[5]).pos[rhop];
                                PP_temp[R][B][Rp][D][rho][beta][rhop][delta]=PP_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPC[iop][A][B][C][D][alpha][beta][gam][delta];


                                // VV
                                for (int mu=1; mu<5; mu++)
                                {
                                Proj_factor=(base_gamma[mu]).entr[rho]*(base_gamma[mu]).entr[rhop];
                                irho=(base_gamma[mu]).pos[rho]; irhop=(base_gamma[mu]).pos[rhop];
                                VV_temp[R][B][Rp][D][rho][beta][rhop][delta]=VV_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPC[iop][A][B][C][D][alpha][beta][gam][delta];
                                }

                                // AA
                                for (int mu=1; mu<5; mu++)
                                {
                                Proj_factor=(base_gamma[mu+11]).entr[rho]*(base_gamma[mu+11]).entr[rhop];
                                irho=(base_gamma[mu+11]).pos[rho]; irhop=(base_gamma[mu+11]).pos[rhop];
                                AA_temp[R][B][Rp][D][rho][beta][rhop][delta]=AA_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPC[iop][A][B][C][D][alpha][beta][gam][delta];
                                }

                                // TT
                                for(int sig=6; sig<12;sig++)
                                {
                                Proj_factor=(base_gamma[sig]).entr[rho]*(base_gamma[sig]).entr[rhop];
                                irho=(base_gamma[sig]).pos[rho]; irhop=(base_gamma[sig]).pos[rhop];
                                TT_temp[R][B][Rp][D][rho][beta][rhop][delta]=TT_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPC[iop][A][B][C][D][alpha][beta][gam][delta];
                                }


                        }// End Internal indices to contract with left propagators

                        SS_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=SS_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + SS_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];
                        PP_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=PP_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + PP_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];
                        VV_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=VV_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + VV_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];
                        AA_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=AA_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + AA_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];
                        TT_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=TT_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + TT_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];


                    } //   End External indices to be contracted with the right propagators

                // add all the intermedial computation for each R Rp rho rhop
                SS[iop]=SS[iop]+SS_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];
                PP[iop]=PP[iop]+PP_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];
                VV[iop]=VV[iop]+VV_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];
                AA[iop]=AA[iop]+AA_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];
                TT[iop]=TT[iop]+TT_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];


                }//End external indices to be contracted with the projector

       /*         printf("iop %d ipm %d\n",iop,ipm);
                printf("SS %8.10f %8.10f \n",creal(SS[iop]),cimag(SS[iop]));
                printf("PP %8.10f %8.10f \n",creal(PP[iop]),cimag(PP[iop]));
                printf("VV %8.10f %8.10f \n",creal(VV[iop]),cimag(VV[iop]));
                printf("AA %8.10f %8.10f \n",creal(AA[iop]),cimag(AA[iop]));
                printf("TT %8.10f %8.10f \n",creal(TT[iop]),cimag(TT[iop]));
       */ 

            } // end iop 

}

void  Projection_TL_OPV( O5_ccccssss OpmPV,colorcolorspinspin Sinv1,colorcolorspinspin Sinv2,colorcolorspinspin Sinv3,colorcolorspinspin Sinv4,int ipm)
{

        ccccssss PS_temp; ccccssss SP_temp; ccccssss AV_temp; ccccssss VA_temp; ccccssss TT5_temp;
        ccccssss PS_temp2; ccccssss SP_temp2; ccccssss AV_temp2; ccccssss VA_temp2; ccccssss TT5_temp2;


	int irho, irhop;
	complex Proj_factor;

            for (int iop=0; iop<5; iop++)
            {
		// Initialize

		SP[iop]=0;
		PS[iop]=0;
		VA[iop]=0;
		AV[iop]=0;
		TT5[iop]=0;
		
		//External indices to be contracted with the projector
		for (int rho=0; rho<4; rho++) for (int rhop=0; rhop<4; rhop++)
		for (int R=0; R<3; R++) for (int Rp=0; Rp<3; Rp++)
		{
                        SP_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;
                        PS_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;
                        VA_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;
                        AV_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;
                        TT5_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=0;

	              // External indices to be contracted with the right propagators	
	              for (int beta=0; beta<4; beta++) for (int delta=0; delta<4; delta++)
	              for (int B=0; B<3; B++) for (int D=0; D<3; D++)
	              {

                        SP_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;
                        PS_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;
                        VA_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;
                        AV_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;
                        TT5_temp[R][B][Rp][D][rho][beta][rhop][delta]=0;

			// Internal indices to contract with left propagators
                        for (int alpha=0; alpha<4; alpha++) for (int gam=0; gam<4; gam++)
                        for (int A=0; A<3; A++) for (int C=0; C<3; C++)
			{
                                // SP
                                Proj_factor=(base_gamma[0]).entr[rho]*(base_gamma[5]).entr[rhop];
                                irho=(base_gamma[0]).pos[rho]; irhop=(base_gamma[5]).pos[rhop]; 
                                SP_temp[R][B][Rp][D][rho][beta][rhop][delta]=SP_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPV[iop][A][B][C][D][alpha][beta][gam][delta];



                                // PS
                                Proj_factor=(base_gamma[5]).entr[rho]*(base_gamma[0]).entr[rhop];
                                irho=(base_gamma[5]).pos[rho]; irhop=(base_gamma[0]).pos[rhop]; 
                                PS_temp[R][B][Rp][D][rho][beta][rhop][delta]=PS_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPV[iop][A][B][C][D][alpha][beta][gam][delta];


                                // VA
                                for (int mu=1; mu<5; mu++)
                                {
                                Proj_factor=(base_gamma[mu]).entr[rho]*(base_gamma[mu+11]).entr[rhop];
                                irho=(base_gamma[mu]).pos[rho]; irhop=(base_gamma[mu+11]).pos[rhop]; 
                                VA_temp[R][B][Rp][D][rho][beta][rhop][delta]=VA_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPV[iop][A][B][C][D][alpha][beta][gam][delta];
				}

                                // AV
                                for (int mu=1; mu<5; mu++)
                                {
                                Proj_factor=(base_gamma[mu+11]).entr[rho]*(base_gamma[mu]).entr[rhop];
                                irho=(base_gamma[mu+11]).pos[rho]; irhop=(base_gamma[mu]).pos[rhop]; 
                                AV_temp[R][B][Rp][D][rho][beta][rhop][delta]=AV_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPV[iop][A][B][C][D][alpha][beta][gam][delta];
				}

                                // TT5
                                for(int sig=6; sig<12;sig++)
                                {
				int sig5=sig+10;
                                Proj_factor=(base_gamma[sig]).entr[rho]*(base_gamma[sig5]).entr[rhop];
                                irho=(base_gamma[sig]).pos[rho]; irhop=(base_gamma[sig5]).pos[rhop]; 
                                TT5_temp[R][B][Rp][D][rho][beta][rhop][delta]=TT5_temp[R][B][Rp][D][rho][beta][rhop][delta]+Proj_factor*Sinv1[R][A][irho][alpha]*Sinv2[Rp][C][irhop][gam]*OpmPV[iop][A][B][C][D][alpha][beta][gam][delta];
				}



			}// End Internal indices to contract with left propagators

                        SP_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=SP_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + SP_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];
                        PS_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=PS_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + PS_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];
                        VA_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=VA_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + VA_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];
                        AV_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=AV_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + AV_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];
                        TT5_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop]=TT5_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop] + TT5_temp[R][B][Rp][D][rho][beta][rhop][delta]*Sinv3[B][R][beta][rho]*Sinv4[D][Rp][delta][rhop];
		    } //   End External indices to be contracted with the right propagators

		// add all the intermedial computation for each R Rp rho rhop

                SP[iop]=SP[iop]+SP_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];
                PS[iop]=PS[iop]+PS_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];
                VA[iop]=VA[iop]+VA_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];
                AV[iop]=AV[iop]+AV_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];
                TT5[iop]=TT5[iop]+TT5_temp2[R][R][Rp][Rp][rho][rho][rhop][rhop];

		}//End external indices to be contracted with the projector

        
	
	    } // end iop 

}


/////////////////////////////////////
// r-prom routines
// not enough memory to do it
/////////////////////////////////////

void Projection_TL_2r_OPC( r_O5_ccccssss OpmPC, colorcolorspinspin *Sinv_x0,colorcolorspinspin *Sinv_0x,int ipm)
{
        //1. Initialize
        for (int iop=0; iop<5; iop++)  SS2r[iop]=PP2r[iop]=VV2r[iop]=AA2r[iop]=TT2r[iop]=0.;

        int Nrcomb=0;
        // 2. Compute and add for each r-comb

                // From eq.26 bibbia
                // Opm(0) = (\psi_1 (x)\bar{\psi}_1(0) ) (\psi_3 (x)\bar{\psi}_3(0) ) Gpm (\psi_2 (0)\bar{\psi}_2(x) ) (\psi_4 (0)\bar{\psi}_4(x) )
                // Opm(0) = S1(x,0) S3(x,0) Gpm S2(0,x) S4(0,x)
                // Lambda = S1(x,0) ^-1 S3(x,0)^-1 Opm S2(0,x)^-1 S4(0,x)^-1

// Como:
// OpmPX(0010)=OpmPX(1011)=OpmPX(1110)*=OpmPX(1000)*
// OpmPX(0001)=OpmPX(1101)=OpmPX(0111)*=OpmPX(0100)*
// (Comprobado numéricamente)
// (Parece lógico porque es cambiar la parte TW - OS pero buscar la simetria) //TODO
// Basta calcular
// OpmPX(0001) y OpmPX(1110) y al promediar quedarnos con la parte real
	  for (int ircomb=0; ircomb<2; ircomb++)
	  {
	  int *rvector=malloc(sizeof(int)*(4)); whichr_sort(ircomb, rvector);
          int r1,r2,r3,r4; r1=rvector[0]; r2=rvector[1]; r3=rvector[2]; r4=rvector[3];
          Projection_TL_OPC(OpmPC[ircomb],Sinv_x0[r1],Sinv_x0[r3],Sinv_0x[r2],Sinv_0x[r4],ipm); 

          for (int iop=0; iop<5; iop++) 
          {
                SS2r[iop]=SS2r[iop]+SS[iop];
                PP2r[iop]=PP2r[iop]+PP[iop];
                VV2r[iop]=VV2r[iop]+VV[iop];
                AA2r[iop]=AA2r[iop]+AA[iop];
                TT2r[iop]=TT2r[iop]+TT[iop];
	

          }
        Nrcomb++;
         }
      // 3. Divide over Nrcomb (with the real part since I'm only computing 0001 & 1110 to save time ) 
      for (int iop=0; iop<5; iop++) 
      {
                SS2r[iop]=creal(SS2r[iop]/Nrcomb);
                PP2r[iop]=creal(PP2r[iop]/Nrcomb);
                VV2r[iop]=creal(VV[iop]/Nrcomb);
                AA2r[iop]=creal(AA2r[iop]/Nrcomb);
                TT2r[iop]=creal(TT2r[iop]/Nrcomb);


    }

}


void Projection_TL_2r_OPV(r_O5_ccccssss OpmPV, colorcolorspinspin *Sinv_x0,colorcolorspinspin *Sinv_0x,int ipm)
{
        //1. Initialize
        for (int iop=0; iop<5; iop++) SP2r[iop]=PS2r[iop]=VA2r[iop]=AV2r[iop]=TT52r[iop]=0.;

        int Nrcomb=0;
        // 2. Compute and add for each r-comb


// Como:
// OpmPX(0010)=OpmPX(1011)=OpmPX(1110)*=OpmPX(1000)*
// OpmPX(0001)=OpmPX(1101)=OpmPX(0111)*=OpmPX(0100)*
// (Comprobado numéricamente)
// (Parece lógico porque es cambiar la parte TW - OS pero buscar la simetria) //TODO
// Basta calcular
// OpmPX(0001) y OpmPX(1110) y al promediar quedarnos con la parte real
          for (int ircomb=0; ircomb<2; ircomb++)
          {
          int *rvector=malloc(sizeof(int)*(4)); whichr_sort(ircomb, rvector);
          int r1,r2,r3,r4; r1=rvector[0]; r2=rvector[1]; r3=rvector[2]; r4=rvector[3];
          Projection_TL_OPV(OpmPV[ircomb],Sinv_x0[r1],Sinv_x0[r3],Sinv_0x[r2],Sinv_0x[r4], ipm);

          for (int iop=0; iop<5; iop++) 
          {

                SP2r[iop]=SP2r[iop]+SP[iop];
                PS2r[iop]=PS2r[iop]+PS[iop];
                VA2r[iop]=VA2r[iop]+VA[iop];
                AV2r[iop]=AV2r[iop]+AV[iop];
                TT52r[iop]=TT52r[iop]+TT5[iop];

          }
        Nrcomb++;
        }
      // 3. Divide over Nrcomb (with the real part since I'm only computing 0001 & 1110 to save time ) 
      for (int iop=0; iop<5; iop++) 
      {


                SP2r[iop]=creal(SP2r[iop]/Nrcomb);
                PS2r[iop]=creal(PS2r[iop]/Nrcomb);
                AV2r[iop]=creal(AV2r[iop]/Nrcomb);
                VA2r[iop]=creal(VA2r[iop]/Nrcomb);
                TT52r[iop]=creal(TT52r[iop]/Nrcomb);
    }

}


void Projection_TL_2r_OpPC(r_O5_ccccssss OpPC,colorcolorspinspin *Sinv_xy, colorcolorspinspin *Sinv_yx)
{
	Projection_TL_2r_OPC(OpPC, Sinv_xy,Sinv_yx,0);
}


